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ABSTRACT 

SIGMA  is  a  set  of  FORTRAN  subprograms  for  solving  the  global  optimization 
problem,  which  implement  a  method  founded  on  the  numerical  solution  of  a 
Cauchy  problem  for  stochastic  differential  equations  inspired  by  quantum 
physics . 

This  paper  gives  a  detailed  description  of  the  method  as  implemented  in 
SIGMA,  and  reports  on  the  numerical  tests  which  have  been  performed  while  the 
SIGMA  package  is  described  in  the  accompanying  Algorithm. 

The  main  conclusions  are  that  SIGMA  performs  very  well  on  several  hard 
test  problems;  unfortunately  given  the  state  of  the  mathematical  software  for 
global  optimization  it  has  not  been  possible  to  make  conclusive  comparisons 
with  other  packages. 
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SIGNIFICANCE  ANT  FXPI..ANA IION 


The  paper  reports  about  a  new  and  very  successful  method  for  finding  a 
"global"  (or  "absolute")  minimum  of  a  functic  i:  of  N  real  variables,  i.e.  the 
point  x  in  N-dimensione.l  none;-  <  ,»>*-•'  ’  ,  ; •.\o  f  the  points)  such  that  not 

only  the  function  j.ncroat>es  one  moves  away  lion  x  in  any  direction, 
i "local"  or  “relative"  mm-rm),  but  disc  such  that  no  other  point  exists 
where  f  has  a  lower  value. 

The  method,  which  was  f ‘ r'p  proposed  by  the  present  authors  in  a  paper 
which  is  to  appear  in  the  Journal  of  Optimization  Theory  and  Applications,  is 
based  on  ideas  from  statistical  mechanics,  and  looks  for  a  point  of  global 
minimum  by  following  the  solution  trajectories  of  a  stochastic  differential 
equation  representing  the  motion  of  a  particle  (in  \-space)  under  the  action 
of  a  potential  field  and  of  a  random  perturbing  force. 

The  paper  gives  a  detailed  description  of  the  complete  algorithm  based  on 
such  a  method,  and  summarize.'-  the  results  of  extensive  numerical  testing  of 
the  FORTRAN  program  implementing  the  algorithm  (the  FORTRAN  program  is 
described  in  a  companion  pap«-r  of  the  same  authors:  Algorithm  SIGMA.  A 
Stochast i c-Integrat i on  Global  Minimization  Algorithm). 

The  tests  have  been  performed  by  running  the  program  on  an  extensive  set 
of  ^aiefuily  selected  test  problems  of  varying  difficulty,  and  the  performance 
ras  been  remarkably  successful,  even  on  very  hard  problems  (e.g.  problems  with 
a  single  point,  of  global  minimum  and  up  to  about  1010  points  of  local 
minimum. ) 


The  method  is  now  being  successfully  tested  on  some  real-world  problems 
in  applied  chemistry,  concerning  the  analysis  of  complex  molecules,  where  one 
iuv.k 6  foi  spatial  patterns  which  are  not  only  stable  (local  minima  of 
potential  energy),  but  have  also  an  absolute  minimum  of  the  potential  energy. 

I  More  generally  there  are  many  problems  in  which  the  solution  depends  on 

the  values  of  several  parameters,  and  the  quality  of  the  solution  can  be 
measured  by  a  single  "performance  figure"  (which  is  therefore  a  function  of 
the  parameters ) ,  e.g.  a  cost,  or  a  loss,  or  a  cost/effectiveness  ratio,  which 
should  be  low,  or  a  gain,  a  utility,  which  should  be  high. 

In  suuh  situations  the  method  can  be  usefully  applied  if  one  is  not 
satisfied  by  finding  a  "sub-ontimal"  solution,  i-e.  a  solution  which  is  the 
best  among  many  other  solutions,  but  one  requires  a  truly  optimal  solution, 
i.e.  the  test  among  all  possible  solutions. 

is  finally  to  be  noted  that,  the  majority  of  the  optimization  methods 
ru  1  >>  ev*1  ilable  deal  with  the  local  optimization  problem,  and  that  no 

methods  of  comparable  power  seem  to  he  available  in  the  field  of  global 

opt.  jin  izatirn. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
;; i -rrm lies  with  MFC,  and  not  with  the  authors  of  this  report. 


A  GLOBAL  OPTIMIZATION  ALGORITHM  USING 
STOCHASTIC  DIFFERENTIAL  EQUATIONS 

Filippo  Aluffi-Pentini^,  Valerio  Parish,  and  Francesco  Zirilli^ 

1.  Introduction. 

In  [1]  a  method  for  solving  the  global  optimization  problem  was 
proposed.  The  method  associates  a  stochastic  differential  equation  with 
the  function  whose  global  minimizer  we  are  looking  for. 

The  stochastic  differential  equation  is  a  stochastic  perturbation 
of  a  "steepest  descent"  ordinary  differential  equation  and  is  inspired 
by  quantum  physics.  In  [1]  the  problem  of  the  numerical  integration  of 
the  stochastic  equations  introduced  was  considered  and  a  suitable  "stochas¬ 
tic"  variation  of  the  Euler  method  was  suggested. 

SICMA  is  a  set  of  FORTRAN  subprograms  implementing  the  above 
method. 

In  sect.  2  we  describe  the  method  as  implemented  in  SIGMA;  in  sect. 

3  we  give  a  general  description  of  the  method  and  some  details  on  the 
implementation;  in  sect.  4  some  numerical  experience  on  test  problems  is 
presented  and  in  sect.  5  conclusions  are  given. 

Unfortunately,  given  the  state  of  the  art  of  mathematical  software 
in  global  optimization,  it  has  not  been  possible  to  make  conclusive  com¬ 
parisons  with  other  packages. 

The  SICMA  package  and  its  usage  are  described  in  the  accompanying 
.Algorithm. 
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2 .  The  method. 

N  N 

Let  R  be  the  N- dimensional  real  euclidean  space  and  let  fdR‘  -*R 

be  a  real  valued  function,  regular  enough  to  justify  the  following  con¬ 
siderations  . 

In  this  paper  we  consider  the  problem  of  'finding  a  global  minimizer 

*  N 

of  £,  that  is,  the  point  x  e  JR'  (or  possible  one  of  the  points)  such 
that 

(2.1)  f(x*)  s  f(x)  Vx  £  IRN 

and  we  propose  a  method  introduced  in  [1]  inspired  by  quantum  physics  to 

compute  numerically  the  global  minimizers  of  f  by  following  the  paths 

of  a  stochastic  differential  equation. 

The  interest  of  the  global  optimization  problem  both  in  mathematics 

and  m  many  applications  is  well  known  and  will  not  be  discussed  here. 

We  want  just  to  remark  here  that  the  root-finding  problem  for  the 

N  N 

system  g(x)  =  0,  where  gdR  -*■  R  can  be  formulated  as  a  global  optimi¬ 
zation  problem  considering  the  function  F(x)  -  ||g(x)i|2>  where  ||  •  ||  2 
is  the  euclidean  norm  in  R^.* 

Despite  its  importance  and  the  efforts  of  many  researchers  the 
global  optimization  problem  is  still  rather  open  and  there  is  a  need  for 
methods  with  solid  mathematical  foundation  and  good  numerical  performance. 


The  present  authors  have  considered  this  idea  both  from  the  mathematical 
point  of  view  (for  a  review  see  [2])  and  from  the  point  of  view  of  producing 
good  software  (see  [3],  [4}).  Tne  method  implemented  in  [3],  [4]  is  in¬ 
spired  by  classical  mechanics,  uses  ordinary  differential  equations,  and  can 
be  regarded  as  a  method  for  global  optimization. 
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Mich  more  satisfactory  is  the  situation  for  the  problem  of  finding 
the  local  minimi zers  of  f ,  where  a  large  body  of  theoretical  and 
numerical  results  exists;  see  for  instance  [5],  [6]  and  the  references 
giver.  therein. 

Or  dinar)’  differential  equations  have  been  used  in  the  study  of  the 
local  optimization  problem  or  of  the  root  finding  problem  by  several 
authors ;  for  a  review  see  [2J. 

The  aho. e  methods  usually  obtain  the  local  minimizers  or  roots  by 
following  the  trajectories  of  suitable  ordinary  differential  equations. 
However,  since  the  property  (2.1)  of  being  a  global  minimizer  is  a  glo¬ 
bal  one.  that  is,  depends  on  the  behaviour  of  f  at  each  point  of  IR^, 
and  the  methods  that  follow  a  trajectory  of  an  ordinary  differential 
equation  are  local,  that  is,  they  depend  only  on  the  behaviour  of  f 
along  the  trajectory,  there  is  no  hope  of  building  a  completely  satis- 
1  sc  eery  method  for  global  optimization  based  on  ordinary  differential 
equations . 

The  situation  is  different  if  we  consider  a  suitable  stochastic 
perturbation  of  an  ordinary  differential  equation  as  explained  in  the 

fol lowing . 

Let  us  first  consider  the  (I to)  stochastic  differential  equation 
•  l.-j  d£  “  -ffCOdt  *  odw 

.on Vf  .  tin  graJk'-ii  cf  f  and  w_('. )  is  a  standard  N- dimensional 

.'.rtr  process,  r.  e  P. . 

aquation  (2.2)  is  known  as  the  Srr.c iuchowski - Kramers  equation  [7]; 
this  equation  i*  a  singular  limit  of  the  Langevin's  equation  when  the 
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The  Smoluchowski -Kramers  equation  has  been  extensively  used  by 
solid  state  physicists  and  chemists  to  study  physical  phenomena  such 
as  atomic  diffusion  in  crystals  or  chemical  reactions. 

In  these  applications  (2.2)  represents  diffusion  across  potential 


barriers  under  the  stochastic  forces  edw,  where  e 


'2kT 

m 


T  is 


the  absolute  temperature,  k  the  Boltzmann  constant,  m  a  suitable  mass 
coefficient,  and  f  is  the  potential  energy. 

We  assume  that 


(2.3)  lim  f(x)  =  +  ® 
li  x  II 2  -*■ 00 

in  such  a  way  that: 

(2.4)  j  e'a2fW  dx  <  »  Ya  «  (]R\{0)) 


and  that  the  minimi zers  of  f  are  isolated  and  non  degenerate. 

It  is  well  known  that  if  j(G(t)  is  the  solution  process  of  (2.2) 
starting  from  an  initial  point  xo  >  the  probability  density  function 
pC(t,x)  of  |_G  (t)  approaches  as  t  -+  »  the  limit  density  pG(x) 
where 


(2.5) 


P„Q0  -  Ae  e 


72  fW 


where  A  is  a  normalization  constant.  The  way  in  which  pE(t,x)  for 
a  class  of  one -dimensional  systems  approaches  pe  (x)  has  been  studied 

oo  — - 

in  detail  by  considering  the  spectrum  of  the  corresponding  Fokker-Planck 
operators  in  [8]  . 
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We  note  that  is  independent  of  x0  anc*  t^iat  as  e  0  P» 
becomes  more  concentrated  at  the  global  minimizers  of  f.  That  is, 

(2.6)  lim  £e(t)  -  5^  in  law 

t-*® 

where  has  a  probability  density  given  by  (2.5)  and 

% 

(2.7)  lim  £ew  -  in  law 

£-+0 

where  is  a  random  variable  having  as  its  probability  density  a 

weighted  sum  of  Dirac's  deltas  concentrated  at  the  global  minimizers  of  f. 
For  example  if  N  *  1  and  f  has  two  global  minimizers  xx>  x2,  with 
(xj)  *  c^  >  0,  i  *  1,2,  we  have  (in  distribution  sense) 

(2.8)  lim  p*(x)  =  y  6(x-xx)  +  (1-y)  5(x-x2) 

e-*0 

where  y  •  (1  +  v/c1/c2  )-1.  In  order  to  obtain  the  global  minimizers  of  f 
as  asymptotic  values  as  t  »  of  a  sample  trajectory  of  a  suitable  sys¬ 
tem  of  stochastic  differential  equations  it  seems  natural  to  try  to  perform 
the  limit  t  -*•  «  (i.e.  (2.6))  and  the  limit  e  -*■  0  (i.e.  (2.7))  together. 

That  is,  we  want  to  consider: 

(2.9)  d£  ■  -Vf (£) dt  +  E(t)dw 
with  i 

(2.10) 

where 


(2.11) 
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In  physical  terms  condition  (2.11)  means  that  the  temperature  T 
is  decreased  to  0  (absolute  zero)  when  t  that  is,  the  system  is 

"frozen". 

Since  we  want  to  end  up  in  a  global  minimizer  cf  f,  that  is,  a 

global  minimizer  of  the  (potential)  energy,  the  system  has  to  be  frozen 

% 

very  slowly  (adiabatically) .  The  way  in  which  e(t)  must  go  to  zero, 
in  order  to  have  that  when  t  the  solution  £(t)  of  (2.9)  becomes 

concentrated  at  the  global  minimizers  of  f,  depends  on  f.  In  par¬ 
ticular,  it  depends  on  the  highest  barrier  in  f  to  be  overcome  to 
reach  the  global  minimizers. 

This  dependence  has  been  studied  using  the  adiabatic  perturbation 
theory  in  [1] .  Similar  ideas  in  the  context  of  combinatorial  optimization 
have  been  introduced  by  Kirkpatrick,  Gelatt,  Vecchi  in  [9], 

In  this  paper  we  restrict  our  attention  to  the  numerical  implemen¬ 
tations  of  the  previous  ideas,  that  is,  the  computation  of  the  global 
minimizers  of  f  by  following  the  paths  defined  by  (2.9),  (2.10).  dis¬ 
regarding  mathematical  problems  such  as  the  difference  between  the  con¬ 
vergence  in  law  of  £(t)  to  a  random  variable  concentrated  at  the  global 
miniiidzers  of  f,  and  the  convergence  with  probability  one  of  the  paths 
of  §_(t)  to  the  global  minimizers  of  f. 

We  consider  the  problem  of  how  to  compute  numerically  these  paths 
keeping  in  mind  that  we  are  not  really  interested  in  the  paths,  but  only 
in  their  asymptotic  values. 

We  discretize  (2.9),  (2.10)  using  the  Euler  method,  that  is  £_(t^) 
i .*  approximated  by  the  solution  of  the  following  finite  difference 


equations : 


/ 


(2.12)  ^  =  +  ^^k^-’k+i  '  -}P  k  =  °'1’2’  ••• 

(2.13)  £0  =  x0 

k- 1 

'-here  t0  =  0 ,  \  *  l  \  >  °,  and  =  wC^) ,  k  =  0,1,2 . 

The  computationally  cheap  Euler  step  seems  a  good  choice  here 
since  in  order  to  obtain  the  global  minimi zers  of  f  as  asymptotic 
values  of  the  paths  e(t)  should  go  to  zero  very  slowly  when  t  ■*  *>, 
and  therefore  a  large  number  of  time  integration  steps  must  be  com¬ 
puted  . 

On  the  right  hand  side  of  (2.12)  we  add  the  random  term 
e(tk)(wk+1  -  wk)  to  the  deterministic  term  -h^  Vf (C^) ,  which  is  com¬ 
putationally  more  expensive  (e.g.  N+l  function  evaluations  if  a 
forward-difference  gradient  is  used),  so  that  the  effort  spent  in  evaluat¬ 
ing  Vf(£.)  is  frequently  lost. 

In  order  to  avoid  this  inconvenience  we  substitute  the  gradient 
? f (£)  with  a  "random  gradient"  as  follows.  Let  r  be  an  N'-dimensional 

random- vector  of  length  1  uniformly  distributed  on  the  N-dimensional  unit 

N 

■sphere.  Then  for  any  given  (non-random)  vector  v  e  1R  its  projection 
along  r  is  such  that: 

(2.11)  N*E(  <r  ,v  >  r)  =  v 


■-here  E(*)  is  the  expected  value,  and  <  •,•  >  is  the  euclidean  inner 
v 

pr.-puct  in  R  . 


Sc  that  in  order  tc  save  numerical  work  (i.e.  functions  evaluations) 
in  (2.12)  we  substitute  7f(£^)  with  the  "random  gradient" 
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(2.15)  Y(£k)  =  N  <  r,  vf(y  >  r. 

1 

i\'e  note  that  since  jt  y($^)  is  the  directional  derivative  in  the  direc¬ 
tion  r,  it  is  computationally  much  cheaper  (e.g.  when  forward  differences 
are  used,  only  2  function  evaluations  are  needed  to  approximate  y(O). 
Therefore,  the  paths  are  computed  approximating*  £(t^)  with  the  solution 
5,  of  the  following  differences  equations: 

4*1  '  5k  "  ‘  V1  k-0.1,2,  ... 

(2.i7)  ;.  ■ 

where  7 (£. )  is  a  finite  difference  (forward  or  central)  approximation 

to  7  (C^) . 

The  complete  algorithm  is  described  in  the  next  section. 
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3.  The  complete  algorithm. 

Ke  give  in  sect.  3.1  a  general  description  of  the  algorithm,  while 
implementation  details  are  given  in  sect.  3.2. 

5.1.  General  description  of  the  algorithm. 

The  basic  time -integration  step  (eq.  (2.16)  and  sect.  3.2.1)  is 
used  to  generate  a  fixed  number  of  trajectories,  which  start  at 

time  zero  from  the  same  initial  conditions  with  possibly  different  values 
of  c(0)  (note  that  even  if  the  starting  values  e(0)  are  equal  the 
trajectories  evolve  differently  due  to  the  stochastic  nature  of  the  inte¬ 
gration  steps) . 

The  trajectories  evolve  (simultaneously  but  independently)  during 
an  "observation  period"  having  a  given  duration  (sect.  3.2.5),  and  with¬ 
in  which  the  noise  coefficient  of  each  trajectory  is  kept  at  a  constant 

value  t  ,  while  the  values  of  the  steplength  h.  and  of  the  spatial 
P  K 

discretization  increment  Ax^  for  computing  the  random  gradient  (eq. 
(2.15)'  and  sect.  3.2.2)  are  automatically  adjusted  for  each  trajectory 
by  the  algorithm  (sects.  3.2.3  and  3.2.4). 

At  the  end  of  every  observation  period  the  corresponding  trajec¬ 
tories  are  compared,  one  of  them  is  discarded  (and  will  not  be  considered 
any  more) ,  all  other  trajectories  are  naturally  continued  in  the  next 
observation  period,  and  one  of  the  trajectories  is  selected  for  "branch¬ 
ing"  (sect.  3.2.6),  that  is  for  generating  also  a  second  continuation 
trajectory  differing  from  the  first  one  only  in  the  starting  values  for 

t  and  Ax,  fsect.  3.2.7),  and  which  is  considered  as  having  the  same 
P  * 

"past  history"  of  the  first  one. 
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Let  Xx  be  the  largest  eigenvalue  of  the  (symmetric  and  non-negative 
definite)  matrix  C. 

We  adopt  the  updating  matrix 
FA  =  6Xx  1  -  C 

where  I  is  the  N*N  identity  matrix,  8  >  1  (B  =  1.3  in  the  present 
implementation),  and  we  obtain  the  updated  valdfe  A*  of  A  by  means  of 
the  formula 

A'  =  aA  F. 

A 

where  a  is  a  normalization  factor  such  that  the  sum  of  the  squares  of 
the  elements  of  A'  is  equal  to  N  (as  in  the  identity  matrix) . 

The  matrix  F^  seems  one  of  the  possible  reasonable  choices, 
since  it  is  positive  definite  for  8  >  1,  it  has  the  same  set  of  eigen¬ 
vectors  as  C,  its  eigenvalue  spectrum  is  obtained  from  the  spectrum  of 

6  \ 

C  by  reflection  around  X  *  ,  and  it  therefore  acts  in  the  right 

direction  to  counter  the  ill -conditioning  of  f. 

The  magnitude  of  the  counter-effect  depends  on  B:  the  adopted 
value  has  been  experimentally  adjusted. 

The  updated  bias  vector  b’  is  chosen  in  order  that  the  scaling 
at  x  does  not  alter  x,  i.e.  in  order  that 
a'x  +  b’  =  Ax  +  b. 

3.2.13  Criteria  for  numerical  equality. 

The  following  two  criteria  are  used  in  a  number  of  places  in  the 
algorithm  to  decide  if  two  given  numbers  x  and  y  are  sufficiently 
close  to  each  other  (within  given  tolerances)  to  be  considered  "numeri¬ 
cally  equal". 
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>.'e  consider  (for  each  trajectory)  the  rescaled  variable  x  =  Ax  +  b, 

.  .■  A  :s  the  rescaling  matrix  and  b  is  a  bias  vector,  and,  instead  cf 
.  a,,  .-.in trice  with  respect  to  x  the  function  f(x)  =  f(x)  =  f(Ax  +  b)  , 

.  *  ry  to  counter  the  ill -conditioning  of  f  with  respect  to  x  by 

;•  ’  ij  justing  A  fursd  b  is  adjusted  in  order  not  to  alter  x)  . 

% 

I.,  updating  cf  A  is  obtained  by  means  of  an  updating  matrix  F^, 

.  If  or  formed  at  the  end  cf  an  observation  period  if  sufficient  data 

..ivu-.uole  jsee  below),  and  if  the  number  of  elapsed  observation  periods 

net  less  than  a  given  number  K  ,  and  greater  than  7N) . 

p3.sca 

Thu  updating  matrix  is  computed  as  described  below,  keeping  in 

...  th.it  the  random  gradients  are  the  only  simply-usable  data  on  the 

u»\  icr  of  f  computed  by  the  algorithm. 

get  >  i  =  1,2,  ...  ,  N,,  be  the  column  vectors  of  the  components 

~  F  C 

.-.I!  .he  X  finite-difference  random  gradients  y  (y  or  y  ) 

:  muted  along  the  trajectory  (also  for  rejected  steps)  from  the  last 


‘ling. 


If  sufficient  data  are  available  (i.e.  if  N  >  2N2)  we  compute  the 


7  =  _L  f  v 

f  N  .7  (i) 

g  i  =  l 

estimated  covariance  matrix 


1  -S  _ .  _  T 


:o  re  a  re. a 


ser.a ble  indicator,  giver,  the  available  data,  of 


1  conditioning  of  f,  as  having  the  larger  eigenvalues 


1 a  ted  with  the  directions  alone  which  the  second  directional  deriva - 


A  e  of  f  is,  on  the  average,  larger. 


. —  ■  '**■  .'V.-i  1.  '  ■  ■  ff»C  ■  V  J.  ...  .  ..... 

21 

We  note  that  each  integration  step  can  be  rejected  only  a  finite 
number  of  times,  each  observation  period  lasts  a  finite  number  of  accepted 
integration  steps,  and  there  is  a  finite  number  of  observation  periods  in 
a  trial;  since  a  finite  number  of  trials  is  allowed,  the  algorithm  will 
stop  after  a  finite  total  number  of  steps  and  of  function  evaluations. 


3.2.11  Admissible  region  for  the  x-values. 

In  order  to  help  the  user  in  trying  to  prevent  computation  failures 
(e.g.  overflow)  the  present  implementation  of  the  method  gives  the  possi¬ 
bility  of  defining  (for  any  given  problem  and  machine  dynamic  range,  and 
based  on  possible  analytical  or  experimental  evidence)  an  admissible  region 
for  the  x-values  (hopefully  containing  the  looked-for  global  minimizer) 
within  which  the  function  values  may  be  safely  computed.  We  use  an 
N- dimensional  interval 


rmin 


<Jti 


1.2. 


.  N, 


where  the  interval  boundaries  must  be  given  before  trial  start. 

Outside  the  admissible  region  the  function  f(x)  is  replaced  by 
an  exponentially  increasing  function,  in  such  a  way  that  the  values  of 
f  and  of  the  external  function  are  matched  at  the  boundary  of  the  region. 


3.2.12  Scaling. 

In  order  to  make  ill-conditioned  problems  more  tractable,  rescaling 
is  performed  by  the  algorithm  as  follows. 


20 


the  preceding  trial,  according  to  the  outcome  (stopping  condition)  of  the 
preceding  trial  and  to  the  number  t  of  trials  performed  from  algorithm 
start,  as  compared  to  the  given  maximum  number  of  trials 
successful  stop:  a  *  103 
unsuccessful  uniform  stop: 

a  =  10  if  t  <  [[(2/5)  Ntrial] ] 
a  =  10'**  otherwise, 

where  [[x]j  is  the  smallest  integer  not  smaller  than  x 
unsuccessful  non-uniform  stop:  a  =  10-4 
The  initial  point  x0  is  selected  as  follows: 
if  t  <[[(2/5)  take  the  value  of  x0  at  algorithm  start 

otherwise  take  Xo  =  *opt 

where  is  the  current  best  minimizer  found  so  far  from  algorithm 

start. 

All  other  initial  values  are  those  of  the  first  trial,  except  the 
initial  values  of  h  and  Aa  which  are  the  values  reached  at  the  end  of 
the  preceding  trial. 


3.2.10  Stopping  criteria  for  the  algorithm. 


The  complete  algorithm  is  stopped,  at  the  end  of  a  trial,  if  a 
given  number  Ng^  has  been  reached  of  uniform  trial  stops  all  at  the 
current  level,  or  in  any  case  if  a  maximum  given  number 

of  trials  has  been  reached. 

Success  is  claimed  by  the  algorithm  if  at  least  one  uniform  stop 
occurred  at  the  current  fQp-p  level. 
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and  the  best  minimum  function  falue  fgpj  found  so  far  from  algorithm 
start:  if  fpj^jN  and  ^Opp  satisfy  at  least  one  of  the  above  criteria, 
with  the  same  tolerances,  the  trial  is  considered  successful  at  the  level 


^OPT’  otherwise  the  trial  is  again  considered  unsuccessful. 

Checking  of  the  stopping  criteria  is  activated  only  if  a  minimum 


given  number  NpMIN 


of  observation  periods  has  been  reached. 


3.2.9  Characteristics  of  the  successive  trials. 

The  operating  conditions  which  are  changed  when  another  trial  is 
started  are: 

-  seed  of  the  random  number  generator 

-  maximum  duration  of  the  trial 

-  policy  for  choosing  tp  for  the  second  continuation  of  a  branched 

trajectory 

-  value  of  e  at  trial  start 

P 

-  initial  point  x0 . 

The  maximum  duration  of  a  trial,  i.e.  the  maximum  number  N~.,v 
of  observation  periods,  is  obtained  as  follows: 

if  the  preceding  trial  had  a  uniform  stop  (sect.  3.2.8)  take  the 
value  of  the  preceding  trial 

otherwise  take  a  value  obtained  by  adding  to  the  preceding  value 
a  fixed  given  increment 

The  policy  for  selecting  ep  for  the  second  continuation  of  a 

branched  trajectory  was  described  in  sect.  3.2.7. 

The  value  of  t  at  the  start  of  a  new  trial  is  obtained  by  means 
P 

of  a  multiplicative  updating  factor  a  applied  to  the  starting  value  of 
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The  updating  factor  for  £p  is  as  follows: 

for  the  first  trial  and  for  any  trial  following  an 


unsuccessful  trial 


F  -  1CT 
e 


where  x  is  a  random  sample  from  a 


standard  normal  distribution 


for  all  other  trials 


F^  =  2y  **  where  y  is  a  random  sample  from  a 

standard  Cauchy  distribution,  i.e.  with  density 

f(y)  =  i/0(i+y2)) 

The  updating  factor  for  Ax^  is: 

3  z 

F,  =10  where  z  is  a  random  sample  from  a  standard 

Ax 

normal  distribution. 


3.2.8  Stopping  criteria  for  a  trial 


A  trial  is  stopped,  at  the  end  of  an  observation  period,  and  after 
having  discarded  the  worst  trajectory,  if  all  the  final  function  values 
of  the  remaining  trajectories  (possibly  at  different  points  x)  are 
"numerically  equal",  i.e.  if  the  maximum,  ^XFH.AX *  an^  t^e 
*TFM1M’  3171011 8  ^ie  final  values  satisfy  at  least  one  of  the  criteria 

in  sect.  3.2.13,  the  relative  difference  criterion  with  a  given  stopping 
tolerance  and/or  the  absolute  difference  criterion  with  given 

stopping  tolerance  ("uniform  stop  at  the  level  fjpMjjyj")  • 

The  trial  is  also  anyway  stopped,  at  the  end  of  the  observation 
period,  if  a  maximum  given  number  •s’p^ay  °f  observation  periods  has  been 
i  cached . 

In  the  latter  case  the  trial  is  considered  unsuccessful,  while  in 


the  former  case  a  comparison  is  made  between  the  final  value  £ 


TFMIN 
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From  the  point  of  view  of  the  noise  coefficient  a  trajectory 
with  larger  ep  is  considered  better  if  the  comparison  is  made  in  an 
early  observation  period  (as  long  as  kp  <  Mp • 1^ ,  where  kp  is  the 
number  of  elapsed  observation  periods,  and  Mp.I^  are  defined  below) 
and  worse  otherwise. 

A  basic  partial  ordering  of  the  trajectories  is  first  obtained  on 
the  basis  of  past  function  values,  and  a  final  total  ordering  is  then  ob¬ 
tained,  if  needed,  by  suitably  exploiting  the  noise-based  ordering. 

The  discarded  trajectory  is  always  the  worst  in  the  ordering,  while 
the  trajectory  selected  for  branching  is  usually  not  the  best  one,  to 
avoid  to  be  stuck  in  a  non-global  minimum. 

Normal  branching  is  performed  on  the  trajectory  which,  in  the  order¬ 
ing,  occupies  the  place  1^  (a  given  integer);  exceptional  branching, 
where  the  best  trajectory  is  selected,  occurs  for  the  first  time  at  the 
end  of  observation  period  kpo,  and  then  every  Mp  periods  (kpQ  and 
Mp  are  given  integers);  i.e.  exceptional  observation  periods  are  those 
numbered 

kp  =kpo+  jMp  Cj  =  °’1’2’  } 

3.2.7  The  second  continuation  of  a  branched  trajectory. 

While  the  first  (unperturbed)  continuation  of  a  trajectory  that 
undergoes  branching  starts  with  the  current  values  of  ep  and  Ax^,  the 
second  continuation  starts  with  values  obtained  by  means  of  multiplica¬ 
tive  random  updating  factors  applied  to  the  current  values. 


In  phase  6a:  y  =  0-1 

We  finally  remark  that  h^  and  Zix^ 
stants  to  avoid  computational  failures. 


are  bounded  by  suitable  con- 


3.2.5  Duration  of  the  observation  period. 

* 

The  duration  of  observation  period  numbered  kp  from  trial  start, 
defined  as  the  number  of  time  integration  steps  in  period  kp,  is 

computed  as  a  function  of  kp  by  means  of  a  formula  which  must  be  chosen 
before  algorithm  start  among  the  following  three  formulas: 


1) 

%  ■  1  *  [l0S2Ckp)) 

("short"  duration) 

2) 

*hp  =  [V 

("medium-size"  duration) 

3) 

Nhp  '  kp 

("long"  duration) 

where 

kp  =  1 , 2 ,  ...  ,  and  [x] 

is  the  largest  integer  not  greater  than 

x. 


3.2.6  Trajectory  selection. 

In  order  to  decide,  at  the  end  of  an  observation  period,  which  tra¬ 
jectory  is  to  be  discarded,  and  which  one  should  be  selected  for  branch¬ 
ing,  we  compare  the  trajectories  on  the  basis  of  the  values  of  their  noise 
coefficient  in  the  observation  period,  and  of  the  function  values  obtained 
from  trial  start. 

From  the  point  of  view  of  past  function  values  a  trajectory  is  con¬ 
sidered  better  than  another  if  it  has  attained  a  lower  function  value  than 
the  other  (excluding  a  possible  initial  part  common  to  both  trajectories) . 
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We  test  fR  and  f^  *  fk  +  Afk  for  numerical  equality  according 
to  the  relative  difference  criterion  (sect.  3.2.13)  with  tolerances 

and  take 


Rl 


=  10'11  and  tr2  =  10"5, 


6  =  2  if  f  and  are  "equal”  within  tR1 

6=)  if  fR  and  fR  are  not  "equal”  within  tr2 

% 

3=1  otherwise. 

The  interval  (10"11,  10" 5)  has  been  adopted  since  it  contains  both 
the  square  root  and  the  cubic  root  of  the  machine  precision  of  most  com¬ 
puters  in  double  precision  (the  square  root  is  appropriate  for  forward 
differences,  while  the  cubic  root  is  appropriate  for  central  differences). 

Updating  factors  y  for  h^ 

In  phase  4a: 

Y  =  1/1.05  for  the  first  attempt  to  the  first  half -step 

Y  =  i  for  the  second  attempt 

Y  =  1/10  for  all  other  attempts 

In  phase  5  the  value  of  y  depends  on  the  current  number  a  of  accepted 
time  integration  steps  in  the  current  observation  period,  and  on  the  cur¬ 
rent  total  number  r  of  half -steps  rejected  so  far  in  the  current  trial 
(excluding  those  possible  rejected  while  attempting  the  first  step). 

If  r  >  0 


Y  =  1 

(if 

a  ^  2r) 

Y  =  1.1 

(if 

2r  <  a  <  3r) 

Y  =  2 

(if 

3r  <  a) 

■  *  0 

Y  =  2 

(if 

a  *  1) 

Y  =  10 

(if 

a  >  1) 
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6a.  If  the  h3lf-step  is  rejected:  reject  also  the  first  half-step, 

update  (decrease)  hk,  and  go  back  to  1. 

6b.  Otherwise:  accept  the  whole  step  and  try  the  next  one. 

Note  however  that  if  the  same  half- step  is  rejected  too  many  times 
the  half-step  is  nevertheless  accepted  in  order*not  to  stop  the  algorithm; 
this  is  not  too  harmful  since  several  trajectories  are  being  computed, 
and  3  "bad"  one  will  be  eventually  discarded  (in  the  present  implementation 
the  bound  is  given  explicitly  for  the  first  half-step  (50  times),  and  im¬ 
plicitly  for  the  second  half-step  (if  1^  becomes  smaller  than  10~so)). 


3.2.4  The  updating  of  h^  and  Ax^. 

The  time -integration  steplength  ^  and  the  spatial  discretization 
increment  Ax^,  for  the  trajectory  under  consideration  are  updated  while 
performing  the  integration  step,  as  described  in  the  preceding  section. 

Updating  is  always  performed  by  means  of  a  multiplicative  updating 
factor  which  is  applied  to  the  old  value  to  obtain  the  new  one. 

The  magnitude  of  the  updating  factors,  as  used  in  the  various 
phases  of  the  sequence  in  the  preceding  sect.  3.2.3,  is  as  follows: 

Updating  factors  8  for  Axk 
In  phase  lb:  S  =  10 6 
In  phase  2a:  8  =10 

In  phase  4b:  6  =  10‘ 4 


In  phase  5  the  value  of  8  depends  on  the  magnitude  of  the  current  esti¬ 
mated  function  increment  Afk  =  |hk|  Ax^  (where  hk  is  hk  or  nk  as 
appropriate),  and  the  function  value  fk  =  f^^. 
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All  attempts  are  with  the  current  (i.e.  updated)  values  of  h^ 
and  Ax^. 

The  sequence  of  attempts  is  as  follows: 

1.  Pick  up  a  random  unit  vector  rk. 

la.  Compute  the  random  increment  s^  (sect.  3.2.2). 

lb.  If  s^  (and  therefore  Ax^)  is  too  small  (i.e.  if  the  square 

of  the  euclidean  norm  of  the  difference  between  the  computed 
values  of  ^  +  and  is  zero,  due  to  the  finite  arithme¬ 
tics  of  the  machine):  update  (increase)  Axk  and  go  back  to  la. 

2.  Compute  (eq.  (3. 2. 2. 2)). 

2a.  If  the  computed  value  of  (\)2  is  zero  (due  to  the  finite 

arithmetics):  update  (increase)  Axk  and  go  back  to  la. 

F 

3.  Compute  the  first  half-step  with  y.  . 

i\ 

Compute  A  'fk  (eq.  (3. 2. 3.1)). 

3a.  If  A'ffc  <  |nkl  Axk 

accept  the  first  half-step  and  jump  to  5. 

r 

4.  Compute  the  first  half -step  with  y^  to  check  the  appropriateness 
Of  toy 

Compute  A'f^  (eq.  (3. 2. 3.1)). 

4a.  If  A'fk  >  Ifif  -  nk|  A^ 

reject  the  half- step,  update  (decrease)  h^  and  go  back  to  1. 

4b.  Otherwise:  accept  the  half- step,  and  update  (decrease)  Ax^. 

5.  Update  (increase)  h^. 

Update  (decrease)  Ax^. 

6.  Compute  the  second  half-step. 

Compute  A,rfk  (eq.  (3. 2. 3. 2)). 


12 


and  the  forward-  and  central -differences  random  gradients 


(3. 2. 2. 3) 


~F  v 
\  =  N  nk 


~C  XT  ~C 

4  =  N  nk  ^ 


F  C 

IVe  use  ^  or  for  y (C^)  in  the  first  half-step  as  described 


in  the  next  section. 


3.2.3  Accepting  and  rejecting  the  half-steps. 


The  computation  of  the  first  half-step  can  be  attempted  with  the 
forward-  or  central -differences  random  gradient  (y^  or  y^  eq.  (3. 2. 2.3)) 
as  described  below. 

In  either  case  the  half-step  is  accepted  or  rejected  according  to 
the  function  increment 


(3. 2. 3.1) 


A'fk  -  f«p  -  fey 


Since  A'fk  should  be  non-positive  for  a  sufficiently  small  value 
of  the  half-step  is  rejected  if  A'fk  is  "numerically  positive", 
i.e.  larger  than  a  given  positive  small  tolerance. 

The  second  half -step  is  rejected  if  the  corresponding  function 
increment 


(3. 2. 5. 2) 


A'\  -  f(^+1)  -  f(£) 


is  positive  and  too  large  (greater  than  100  e*  in  the  present  implemen¬ 
tation)  . 

The  sequence  of  attempts  affects  the  updating  of  h^  and  Ax^  as 
described  below;  the  amount  of  the  updating  is  described  in  sect.  3.2.4. 
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The  basic  step  (3. 2. 1.1)  is  actus”''  performed  in  two  half -steps 

(3. 2. 1.2)  ^  *  £k  -  hk  £0^.)  (first  half-step) 

and 

(3. 2. 1.3)  C^+1  =  ^  *  £p^Hk  -  (second  half-step) 

Both  half -steps  depend  on  hk  while  the  first  depends  also  on  the 
current  value  Axk  of  the  spatial  discretization  increment  used  in  com¬ 
puting  y(Ck). 

Either  half-step  can  be  rejected  if  deemed  not  satisfactory,  as 
described  in  sect.  3.2.3. 


3.2.2  The  finite-differences  random  gradient. 

Given  the  current  value  Axk  of  the  spatial  discretization  incre¬ 
ment  for  the  trajectory  under  consideration,  we  consider  the  random  in¬ 
crement  vector 


4  =  ^'4 


where  r,  is  a  random  sample  of  a  vector  uniformly  distributed  on  the 
unit  sphere  in  R  ,  the  forward  and  central  differences 


(3.2. 2.1) 


aFfk  -  «4  *  4)  • 

A  4  *  41  ’  "  4)J 


the  forward-  and  central-differences  directional  derivatives 
(5. 2. 2. 2)  nk  =  AFfk/Axk 


The  set  of  simultaneous  trajectories  is  considered  as  a  single  trial, 
which  is  stopped  as  described  in  sect.  3.2.8,  and  is  repeated  a  number  of 
times  with  different  operating  conditions  (sect.  3.2.9). 

The  stopping  criteria  for  the  complete  algorithm  are  described  in 
sect.  3.2.10. 

% 

The  use  of  an  admissible  region  for  the  x-values  is  described  in 
sect.  3.2.11,  scaling  is  described  in  sect.  3.2.12,  and  criteria  for 
numerical  equality  in  sect.  3.2.13. 

3.2.  Implementation  details. 

3.2.1  The  time -integration  step. 

The  basic  time -integration  step  (eq.  (2.16))  is  used,  for  the  tra¬ 
jectory'  under  consideration,  in  the  form 

(3.2. 1.1)  4*1  '  k  -  \  i'V  +  ep  \  \  O'  -  ••• ) 

where  h,  and  e  are  the  current  values  of  the  steplength  and  of  the 

K  p 

noise  coefficient  (the  noise  coefficient  has  a  constant  value  through¬ 
out  the  current  observation  period  (sect  3.1));  u^  is  a  random  vector 
sample  from  an  N-dimensional  standard  Gaussian  distribution,  and 

•S.  "k  *  Vi 

due  to  the  properties  of  the  Wiener  process. 

The  computation  of  the  finite-differences  random  gradient 
is  described  in  the  next  section. 


a)  Relative  difference  criterion 


|x-y |  <  (|x|  +  |y|)/2 

b)  Absolute  difference  criterion 

|x-y!<tABs 
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4.  Numerical  Testing. 

SIGMA  has  been  numerically  tested  on  a  number  of  test  rpoblems  run 
on  two  computers.  The  test  problems  are  described  in  sect.  4.1,  the  com¬ 
puters  in  sect.  4.2  and  some  numerical  results  arereported  in  sect.  4.3. 

% 

4.1.  Test  problems. 

The  set  of  test  problems  is  fully  described  in  [19]  together  with 
the  initial  points;  the  test  problems  are; 

1.  A  fourth  order  polynomial  (N  *  1) 

2.  Goldstein  sixth  order  polynomial  (N  *  l.) 

3.  One  dimensional  penalized  Shubert  function  (N  =  1] 

4.  A  fourth  order  polynomial  in  two  variables  (N  «*  2) 

5.  A  function  with  a  single  row  of  local  minima  (N  *  2) 

6.  Six  hump  camel  function  (N  =  2) 

7.  Two  dimensional  penalized  Shubert  function  $  =  0  (N  *  2) 

8.  Two  dimensional  penalized  Shubert  function  6  =  0.5  (N  =  2) 

9.  Two  dimensional  penalized  Shubert  function  6=1  (N  =  2) 

10.  A  function  with  three  ill-conditioned  minima  a  *  10  (N  =  2) 

11.  A  function  with  three  ill-conditioned  minima  a  =  100  (N  =  2) 

12.  A  function  with  three  ill-conditioned  minima  a  *  1000  (N  =  2) 

13.  A  function  with  three  ill-conditioned  minima  a  =  10000  (N  =  2) 

14.  A  function  with  three  ill-conditioned  minima  a  =  10s  O'  =  2) 

15.  A  function  with  three  ill-conditioned  minima  a  =  106  (N  =  2) 

16.  Goldstein-Price  function  O’  =  2) 

17.  Penalized  Branin  function  (N  =  2) 

18.  Penalized  Shekel  function  M  =  5  (N  =  4) 
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19.  Penalized  Shekel  function  M  =  7  CN  =4) 

20.  Penalized  Shekel  function  M  =10  (N  =  4) 

21.  Penalized  three  dimensional  Hartman  function  (N  =  3) 

22.  Penalized  six  dimensional  Hartman  function  (N  =  6) 

2.3.  Penalized  Levy  Montalvo  function,  type  1  (N  =  2) 

% 

24.  Penalized  Levy  Montalvo  function,  type  1  (N  =  3) 

25.  Penalized  Levy  Montalvo  function,  type  1  (N  =  4) 

26.  Penalized  Levy  Montalvo  function,  type  2  (N  =  5) 

27.  Penalized  Levy  Montalvo  function,  type  2  (N  =  8) 

28.  Penalized  Levy  Montalvo  function,  type  2,  (N1  =  10) 

29.  Penalized  Levy  Montalvo  function,  type  3,  range  10  (N  =  2) 

30.  Penalized  Levy  Montalvo  function,  type  3,  range  10  (N  *  3) 

51.  Penalized  Levy  Montalvo  function,  type  3,  range  10  CN  *  4) 

32.  Penalized  Levy  Montalvo  function,  type  3,  range  5  (N  =  5) 

33.  Penalized  Levy  Montalvo  function,  type  3,  range  5  (N  =  6) 

34.  Penalized  Levy  Montalvo  function,  type  3,  range  5  (N  =  7) 

35.  A  function  with  a  cusp  shaped  minima  (N  =  5) 

36.  A  function  with  a  global  minimum  having  a  small  region 

of  attraction  a  =  100  (N  =  2) 

37.  A  function  with  a  global  minimum  having  a  small  region 

of  attraction  a  =  10  (N  =  5) 

We  used  the  above  functions,  and  the  standard  initial  points  as 
i re  coded  in  the  subroutines  GLCMTF  and  GLQMIP,  which  are  available 
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4.2. '  Test  computers. 

We  considered  two  typical  machines  of  "large"  and  "small"  dynamic 
range,  that  is,  with  11  and  8  bits  for  the  exponent  (biased  or  signed) 
of  double  precision  numbers,  and  corresponding  dynamic  range  of  about 

4*  3  0  6  +38 

10“  and  10“  .  The  tests  were  actually  performed  on: 

♦ 

—  UNIVAC  1100/82  w'ith  EXEC8  operating  system  and  FORTRAN  (ASCII) 
conputer  (level  10R1)  ("large"  dynamic  range) 

—  D.E.C.  VAX  11/750  with  VMS  operating  system  (vers.  3.0) 
and  FORTRAN  compiler  (vers.  3)  ("small"  dynamic  range) 

4.3.  Numerical  results. 

Numerical  results  of  running  SIGMA  on  the  above  problems  and  on  the 
above  machines  are  described  below.  All  results  were  obtained  under  the 
following  operating  conditions. 

The  easy-to-use  driver  subroutine  SIGMA1  (described  in  the  accompany¬ 
ing  algorithm)  was  used,  with  *  1,2, 3, 4, 5.  All  numerical  values  used 

for  the  parameters  are  set  in  the  driver  SIGMA1  and  in  the  other  subroutines 
which  are  described  in  the  accompanying  Algorithm. 

All  numerical  results  are  reported  on  Tables  1,  2,  and  3.  Table  1 
reports  some  performance  data  (i.e.  output  indicator  IOUT  and  number  of 
functions  evaluations)  as  obtained  from  SIG^A  output  for  each  of  the  37 
test  problems  and  for  the  testing  both  on  the  "large"  and  "small"  dynamic 
range  machines.  In  order  to  evaluate  the  performance  of  SIGMA  we  consider 
all  the  cases  in  which  the  program  claimed  a  success  (output  indicator 
IOUT  >0)  or  a  failure  (IOUT  <  0)  and  —  by  comparing  the  final  point 
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with  the  known  solutions  —  we  identify  the  cases  in  which  such  a  claim 
is  clearly  incorrect  (i.e.  success  claim  when  the  final  point  is  not  even 
approximately  close  to  the  known  solution,  or  failure  claim  when  the  final 
point  is  practically  coincident  with  the  known  solution) .  It  is  also 

meaningful  to  consider  all  the  cases  in  which  a  computational  failure 

% 

due  to  overflow  actually  occurrs  at  any  point  of  the  iteration. 

Table  2  and  Table  3  report  for  each  problem  and  summarized  for  all 
problems  data  concerning  the  effectiveness ,  dependability  and  robustness 
—  in  the  form  of  total  numbers  of  correctly  claimed  successes,  correctly 
claimed  failures,  incorrect  success  or  failure  claims  and  total  number  of 
overflows  —  for  the  two  machines. 


.  -VlS-  i  _• _ .  .  ' 


UN  I  VAC 


29 


o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

© 

o 

o 

<*~s 

o 

© 

o 

o 

© 

o 

O 

rH 

LO 

rH 

to 

rH 

o 

LD 

SO 

rH 

LTi 

cc 

nC 

LO 

r^ 

CM 

© 

o 

o 

rN 

CM 

(NJ 

LO 

rH 

CM 

to 

CD 

CD 

sO 

LO 

LO 

LO 

rH 

o 

LO 

© 

L0 

CM 

lO 

CM 

© 

00 

o 

CD 

in 

o 

o 

© 

to 

o 

© 

o 

oo 

CD 

CM 

LO 

o 

to 

to 

CC 

cr. 

to 

to 

to 

rH 

oc 

N4 

c- 

rH 

00 

CD 

(NJ 

rH 

sO 

sO 

r- 

to 

n* 

00 

LO 

© 

CM 

o 

C- 

© 

CM 

CM 

LO 

N. 

© 

o 

© 

in 

to 

CD 

o 

sO 

LO 

HJ- 

to 

rH 

00 

cc 

o 

(NJ 

00 

LO 

cr. 

rH 

rH 

rH 

to 

in 

% 

rH 

(NJ 

to 

to 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

© 

o 

© 

© 

© 

o 

© 

o 

© 

o 

o 

o 

o 

o 

c- 

CD 

to 

00 

to 

(SI 

in 

to 

00 

LO 

00 

CM 

c- 

LO 

LO 

rN 

oo 

LO 

cc 

NJ 

rH 

oo 

OO 

CM 

Nl 

co 

sO 

lO 

o 

rH 

to 

CM 

oc 

© 

IN 

rN 

CM 

oc 

00 

o 

LO 

© 

rH 

in 

CD 

SO 

in 

rN 

oo 

to 

00 

to 

LO 

CC 

(NJ 

rH 

© 

to 

© 

© 

rH 

CM 

CM 

in 

CD 

lO 

tN 

CD 

CD 

CD 

rH 

© 

© 

[n 

rH 

to 

rN 

CC 

CT) 

TT 

oo 

o 

to 

to 

LO 

© 

CD 

LO 

lO 

LO 

CD 

sO 

(NJ 

CM 

to 

CM 

o 

CM 

o 

(NJ 

IN 

© 

to 

OO 

to 

sO 

CM 

NJ 

CM 

o 

o 

© 

O 

O 

o 

o 

o 

O 

o 

© 

o 

o 

o 

o 

o 

O 

O 

© 

o 

o 

o 

o 

n- 

to 

CD 

CM 

o 

to 

rH 

o 

CO 

© 

(NJ 

lO 

CM 

in 

rH 

00 

co 

© 

© 

CM 

00 

cr> 

sO 

CD 

rH 

00 

© 

© 

LO 

o 

to 

LO 

© 

© 

o 

© 

CM 

© 

oo 

CM 

LO 

rH 

IN 

IN 

O 

CO 

00 

lO 

to 

so 

rH 

to 

rH 

CD 

rH 

© 

© 

rH 

(NJ 

to 

as 

to 

© 

*3* 

to 

LO 

fO 

O 

to 

o 

oo 

CD 

CD 

CO 

o 

CN) 

to 

00 

© 

LO 

© 

00 

rH 

rH 

CM 

rN 

o 

NJ 

CM 

CM 

nj 

K1 

to 

rH 

in 

rH 

(NJ 

rH 

(NJ 

rH 

rH 

© 

< — 4 

o 

to 

© 

NJ 

CM 

to 

to 

rH 

rH 

NJ 

© 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

O 

o 

o 

© 

o 

o 

o 

O 

o 

in 

CD 

rH 

00 

rH 

to 

© 

CD 

LO 

rH 

© 

[N 

rH 

o 

o 

o 

© 

CM 

o 

CO 

o 

o 

o 

cn 

NJ 

rf 

LO 

to 

CM 

to 

OO 

© 

00 

LO 

to 

© 

IN 

LO 

LO 

cn 

NJ 

co 

LO 

tN 

oo 

CM 

to 

CNJ 

o 

© 

In 

to 

(NJ 

(NJ 

© 

O 

lO 

© 

LO 

LO 

o 

LO 

NT 

o 

rH 

CD 

in 

LO 

to 

00 

LO 

in 

to 

O 

cc 

in 

o 

o 

rH 

o 

© 

o 

^3* 

o 

NJ 

O 

o 

rH 

rH 

rH 

cm 

to 

to 

oo 

rH 

rH 

rH 

(NJ 

-H 

LT) 

© 

rH 

rH 

rH 

o 

o 

O 

c 

o 

o 

o 

o 

o 

O 

o 

o 

O 

© 

o 

o 

© 

© 

© 

o 

o 

o 

O 

oo 

oo 

o 

in 

LO 

o 

© 

sO 

rH 

© 

LO 

rH 

r- 

rH 

o 

o> 

LO 

o 

c*> 

00 

oo 

LO 

to 

oo 

CT> 

OO 

o 

00 

cr> 

Cl 

© 

o 

tN 

(NJ 

© 

rN 

LO 

LO 

CM 

LD 

CM 

© 

to 

sO 

sD 

rH 

SO 

rH 

LO 

rH 

n- 

rH 

© 

o 

(N 

o 

tN 

CTj 

*3* 

IN 

oo 

tO 

to 

00 

o 

nj 

04 

CM 

LO 

SO 

to 

to 

in 

..  r» 

in 

o 

CM 

CT) 

to 

rH 

to 

rH 

rH 

rH 

rH 

rH 

NI 

rsj 

NJ 

cm 

CM 

(NJ 

CM 

(NJ 

CM 

CM 

r« 

NJ 

CNJ 

(NJ 

rr 

rr 

to 

NJ 

i— < 

N] 

to 

lO 

o 

in 

00 

CD 

o 

*H 

NJ 

rO 

lO 

fN 

oc 

o 

i-H 

NJ 

K'J 

rH 

*H 

'  ‘ 

rH 

r_— 

r—i 

r-^ 

rH 

M 

NJ 

NJ 

NJ 

NJ 

k  ■*  -  '  V  "  "  >  ”.*•  "  ■  "  *  "'^*T  ’  w-  V*  1  *T  •  J  •  t  w  1 


30 


P 

B 


,  -— 

» 

j 

O 

o 

o 

o 

o 

o 

O 

o 

o 

o 

o 

o 

o 

eg 

r— 1 

r- 

cc 

ro 

to 

eg 

to 

r— ( 

re 

rH 

i  ^ 

eg 

i — i 

tO 

rH 

rg 

r- 

vC 

re 

rr 

re 

\D 

CO 

o 

i :  — 

,  — > 

r- 

o 

CC 

rg 

o 

<3- 

lO 

o 

OO 

re 

re 

i  X- . 

i  ~ 

- 

re 

rr 

<c 

r- 

O 

o 

C 

o 

to 

ci 

rg 

i 

i 

c- 

G 

lO 

o 

oo 

o 

r*- 

l n 

rg 

CO 

cn 

: 

— * 

re 

O 

rr 

r— i 

rg 

rg 

o 

re 

LO 

rH 

\ 

t 

1 

j-'X  I  c—‘ 

-• 

O 

o 

o 

© 

o 

o 

o 

O 

O 

O 

o 

o 

i 

1 

i  -. 

03 

'~S 

fO 

r- 

re 

o 

Cl 

rH 

rH 

CTi 

rg 

o 

O' , 

rg 

re 

o 

to 

oo 

o 

•^r 

Cl 

to 

to 

•-►  1 

Cl 

r— < 

Ci 

rO 

re 

© 

oo 

o 

rg 

CTj 

o 

sO 

z 

K* 

•> 

r 

•» 

*■ 

•* 

•* 

•• 

A 

r 

• 

I  ° 

r- 

O’ 

LO 

re 

pmmd 

r- 

re 

r- 

to 

Ci 

rg 

1 

lO 

to 

r- 

r-H 

lO 

o 

iO 

o 

o 

rg 

r- 

o 

j 

rr 

vO 

re 

*-H 

rg 

rg 

ro 

ro 

Tj- 

r-H 

1 

< 

o  i  ... 

i  w 

o 

~ 

o 

o 

o 

O 

O 

O 

o 

o 

O 

o 

o 

'-i 

i 

un 

o 

cc 

CO 

rr 

o 

rg 

rg 

00 

lO 

to 

LO 

— >  i 

• — < 

LO 

-n 

rg 

o 

'rr 

*-h 

LO 

o 

oo 

r- 

ro 

.  c 

G 

o 

r- 

OO 

cr> 

Cl 

rg 

CO 

o 

oo 

rH 

rg 

rH 

_ _ 

^  , 

•* 

•" 

•* 

* 

* 

*• 

*■ 

•* 

.  r  - 

r— ' 

wO 

rr 

CC 

r- 

CO 

LO 

rg 

r- 

o 

to 

LO 

r-*- 

r  « 

■C-f 

r- 

LO 

cr. 

o 

r- 

r^» 

LO 

r- 

rg 

rg 

• 

1 

to 

re 

rg 

rg 

rg 

ro 

►-•  1 

1 

- 

~ 

- 

o 

© 

- 

O 

O 

O 

O 

o 

o 

O 

1  i- 

r  j 

rr 

f-H 

LO 

cr. 

r- 

o 

rg 

o 

00 

o 

Ci 

! 

uO 

■ — ■ 

O 

r- 

r- 

cc 

lO 

LO 

LO 

O 

i* 

i 

r  g 

rH 

•-H 

sQ 

r^ 

r- 

oo 

o 

rg 

to 

rH 

i  - 

\  -  ^ 

~T 

rH 

o 

LO 

ro 

ro 

o 

rg 

00 

o 

o 

rg 

1  m 

L-O 

cc 

rg 

rg 

to 

f-H 

o 

Ci 

LO 

rH 

rH 

1 

r  i 

rr 

f-H 

rH 

rH 

rH 

rg 

i 


l 


G  ' 

° 

o 

G 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

! 

!c. 

to 

rr 

LO 

rg 

re 

o 

r-H 

oo 

LO 

rg 

00 

*— *  i  c» 

LO 

00 

rg 

r- 

o 

to 

rH 

rH 

m 

to 

G 

ro 

^ 

Ci 

o 

to 

O 

re 

r- 

LO 

Ci 

f-H 

LO 

o 

rH 

o 

!  « 

*-H 

to 

^4 

o 

LO 

LO 

rH 

O 

to 

rj* 

rO 

ro 

vO 

*~1 

re 

LTi 

~g 

r-^ 

•— 1 

LO 

LO 

rH 

to 

7!  I  i^.  re  l.o  oo  o  eg  to  re  tn  o  r^  i/i  m  l/i 


^  j  ^  i 

lO  o  M  0*1  o  H  N  to  re  1/1  o  r» 

/.  I  o.  In  rj  N  N  rg  rg  ro  ro  to  rO  tO  rO  ro  tO 

1  ^  i 


r 


Table  1  (continued) 


h- 


r 


N 


y 

to 

tO 

Oi 

tO 

a> 

OO 

rH 

(NJ 

(NJ 

lO 

tO 

in 

rH 

o 

rH 

o 

l o 

CT> 

rs. 

(NJ 

o 

Ct 

lo 

LO 

csi 

«— < 

ro 

a i 

o 

o 

y 

rH 

at 

to 

rs) 

to 

c 

c 

r^. 

o 

y 

Oi 

Ct 

y 

to 

o 

rsj 

y 

00 

LO 

o 

00 

r- 

to 

y 

tH 

oo 

to 

y 

cs 

to 

vO 

o 

o 

rs) 

o 

y 

l n 

rH 

(NJ 

at 

at 

CNJ 

rH 

r— • 

o 

to 

LO 

o 

oo 

LO 

o 

o 

tO 

lo 

o 

Csi 

Ct 

y 

y 

to 

to 

*3- 

o 

to 

cr. 

(NJ 

LO 

to 

rK 

r— ( 

rH 

to 

(NJ 

rH 

<nj 

rsj 

tH 

T-H 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

© 

o 

© 

o 

o 

© 

to 

(NJ 

r-H 

y 

r- 

at 

y 

to 

at 

o 

oo 

at 

at 

CO 

r- 

LO 

LO 

tN» 

o 

o 

(SJ 

OO 

cc 

(NJ 

y 

tO 

rH 

00 

LO 

O) 

o 

LO 

00 

(NJ 

y 

to 

o 

(NJ 

y 

CSI 

y 

00 

oo 

r-. 

o 

to 

rH 

(NJ 

o 

00 

y 

to 

y 

y 

r- 

to 

lO 

to 

r-H 

to 

at 

CO 

to 

o 

LO 

r-H 

r- 

r- 

o 

at 

o 

rH 

LO 

(NJ 

oo 

00 

(NJ 

to 

o 

o 

rH 

C7t 

at 

oo 

CSJ 

y 

LO 

to 

(NJ 

lo 

CO 

LO 

LO 

to 

CNJ 

to 

y 

(NJ 

(NJ 

(NJ 

(NJ 

tO 

oc 

ift 

tO 

o 

fN 

OO 

L0 

tO 

lH  Lf)  iH  N  H 


<0 


o  o 


o  o 


c  o 


o  o  o 


to 

rH 

r- 

CN 

to 

CSJ 

to 

LO 

tO 

fO 

LO 

rH 

00 

at 

rH 

LO 

y 

y 

oc 

r>- 

to 

y 

(SJ 

-a- 

fN- 

LO 

© 

(NJ 

tN. 

tO 

at 

at 

tO 

at 

at 

o 

rH 

to 

LO 

y 

o 

at 

-r 

Ct 

y 

to 

rH 

to 

y 

oo 

rH 

LO 

to 

tN. 

rH 

LO 

r>- 

to 

o 

y 

(NJ 

rH 

a. 

oo 

LO 

o 

r-H 

-a- 

o 

OC 

© 

LO 

(SJ 

at 

rH 

at 

(SJ 

00 

LO 

t*^ 

y 

oc 

at 

y 

-T 

to 

(SJ 

(NJ 

rH 

(NJ 

(NJ 

to 

CNJ 

00 

00 

r*- 

rH 

rH 

rH 

r-l 

rH 

rH 

(SJ 

rH 

(NJ 

y 

(SJ 

y 

rH 

LO 

<u 


oooooooooooooooo 


o  o  o  o 


► 

t 

c 

► 


Jk 


f-. 

(SI 

(SJ 

to 

tO 

o 

00 

00 

LO 

at 

rH 

at 

rH 

to 

rH 

rH 

o 

(NJ 

rH 

(NJ 

tO 

to 

LO 

y 

y 

rH 

rH 

y 

o 

y 

CO 

LO 

o 

y 

to 

HJ* 

00 

at 

o 

LO 

LO 

*-H 

rH 

rH 

CO 

tO 

.  LO 

tO 

C^- 

o 

y 

to 

co 

a*. 

at 

Ct 

O 

y 

to 

y 

o 

(NJ 

at 

to 

OO 

tO 

to 

(NJ 

to 

LO 

r- 

at 

rH 

y 

to 

LO 

at 

(N- 

at 

rH 

at 

to 

o 

rH 

© 

or 

rH 

00 

a> 

rH 

rH 

rH 

rH 

rH 

y 

to 

r-H 

rH 

rH 

rH 

to 

rH 

rH 

CNJ 

(SJ 

o 


II 


'sJ 

to 

"y* 


«*r» 

9 

a. 

z 


° 

o 

o 

o 

o 

o 

o 

O 

o 

o 

o 

o 

o 

o 

o 

o 

c 

o 

o 

o 

o 

o 

o 

rsj 

rH 

tO 

LO 

y 

o 

to 

00 

rH 

^H 

00 

to 

CT> 

o 

to 

OO 

rH 

y 

LO 

y 

(NJ 

to 

(NJ 

tO 

at 

LO 

y 

y 

rH 

LO 

rH 

O 

rH 

to 

a 

r- 

y 

to 

IN. 

to 

LO 

LO 

(NJ 

o 

o 

to 

o 

to 

r-H 

rsj 

rsi 

to 

y 

^y 

r^- 

to 

at 

LO 

r-- 

to 

-H 

a. 

(NJ 

y 

o 

y 

•«y 

to 

o 

tO 

o 

tO 

lO 

(NJ 

(NJ 

•y 

y 

y 

to 

LO 

rH 

rH 

rH 

(NJ 

(N| 

rH 

rH 

r-H 

(NJ 

(NJ 

(NJ 

rsj 

(NJ 

(SJ 

(NJ 

(SJ 

(SJ 

(NJ 

(NJ 

(NJ 

(NJ 

(NJ 

y 

y 

to 

o> 

(NJ 

rH 

(SJ 

to 

Tf 

LO 

tO 

oo 

at 

o 

^H 

(SI 

to 

y 

lO 

o 

f'N 

00 

at 

rH 

(NJ 

to 

rH 

»—H 

rH 

rH 

rH 

rH 

-H 

*— H 

(NJ 

(NJ 

(NJ 

(NJ 

.1 


VAX  (continued) 


32 


oo 

s 

c 

•H 

•  rH 

*3 

•3 

3 

o 

o 

o 

o 

o 

o 

O 

o 

o 

o 

o 

© 

o 

O 

3 

B 

U 

c 

•  H 

e 

o 

o 

CN) 

Cl 

CN) 

CN) 

LO 

LO 

sO 

CN) 

o 

Ci 

o 

00 

l/> 

"3 

rvi 

OO 

Cl 

CN) 

cn 

r- 

© 

o 

o 

00 

»H 

vO 

c 

C 

Cl 

o 

CN) 

rH 

Ci 

rH 

o 

OO 

Cl 

LO 

CN) 

o 

3 

•> 

** 

* 

fh 

A 

A 

•v 

* 

a 

•* 

0S 

•  rH 

i-t 

o 

UO 

o 

rH 

fO 

vO 

1 0 

C C 

© 

CN) 

c- 

4-> 

• 

cn 

CNJ 

oo 

CO 

rN. 

to 

Cl 

rH 

Os 

■*r 

o 

OO 

o 

C8 

rH 

Cs) 

in 

ro 

•*T 

rH 

to 

to 

to 

CN) 

to 

Cs) 

*“"*  * 

3 

o 

rH 

j= 

rt 

H 

> 

o 

o 

4-1 

o 

o 

o 

o 

o 

o 

o 

o 

O 

o 

o 

o 

o 

o 

c 

3 

o 

p. 

•  rH 

4-> 

o 

" 

U 

u 

o 

00 

UO 

oo 

LO 

00 

LO 

LO 

rH 

rH 

o 

c-n 

cn 

O 

3 

o 

Cn) 

01 

rv) 

rH 

(71 

rH 

© 

to 

C- 

00 

00 

oo 

cn 

rH 

HH 

4-1 

LTi 

CO 

o 

Cl 

LO 

to 

o 

o 

Cl 

in 

• 

•• 

0s 

•t 

0s 

0s 

•* 

* 

0s 

* 

0s 

0S 

* 

Cs) 

4h 

-3 

cn 

Ci 

vO 

o 

rH 

rH 

Cn) 

in 

Cl 

m 

o 

• 

o 

0) 

o 

CN) 

T - < 

l n 

fO 

o 

oo 

oo 

o 

Cn) 

m 

no 

'O 

rH 

tj- 

fO 

rH 

CN) 

CN) 

CN) 

rH 

Cn) 

4. 

3 

• 

u 

3 

u 

>9 

S 

o 

3 

V) 

tfl 

c 

3 

c 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

e 

o 

rH 

o 

o 

rt 

</> 

4-1 

3 

II 

o 

4-> 

5 

C- 

LO 

Ol 

CN) 

rr 

c- 

to 

o 

!"■* 

in 

in 

m 

rH 

o 

OO 

c- 

cn 

LO 

in 

CN) 

CO 

© 

00 

i/i 

to 

LO 

rH 

2 

II 

in 

(NJ 

o 

00 

Cn) 

Ci 

rH 

OO 

Ol 

cn 

(M 

4 -1 

A 

0S 

«s 

#v 

* 

« 

r 

A 

0* 

0t 

0> 

£ 

O 

o 

© 

CN) 

rH 

LO 

00 

00 

© 

Cl 

m 

Cl 

to 

Cl 

Z 

OO 

LD 

© 

o 

rH 

cn 

rH 

00 

m 

t'' 

CN) 

(S 

to 

rH 

rH 

rH 

rH 

t-H 

36 


Table  2  (continued) 

VAX  (continued) 

NSUC  _ _1 _ 2 _ 3 _ 4 _ 5_ 

XPROB  N 

29  2  1  1  1  1  1 

30  3  1  1  1  1  1 

31  4  1  ‘  1  1  1  1 

32  5  1  1  1  1  1 

33  6  1  1  1  1  1 

34  7  1  1  1  1  1 

35  5  1  1  1  '  1  1 

36  2  3  3  3  3  3 

37  5  3  3  3  3  3 

1  =  success  correctly  claimed 

2  =  failure  correctly  claimed 

3  =  incorrect  claim 


4  =  overflow 
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Conclusions . 


The  SIGMA  package  presented  here  seems  to  perform  quite  well  on 

too  proposed  test  problems. 

As  it  is  shown  in  [10]  some  of  the  test  problems  are  very  hard; 
r  example,  Problem  28  (N  =  10)  has  a  single  global  minimi zer  and  a 
cu'.vrr  v£  local  minimi zers  of  order  101C  in  the  region  |x^|  <  10 

i  »  1,2,  ....  10. 

i'abie  2  shows  that  from  the  point  of  view  of  the  effectiveness  as 
red  by  the  number  of  correctly  claimed  successes  the  performance 

of  liGMA  is  very  satisfactory;  moreover,  it  is  remarkably  machine  inde- 

# 

:u'V.;ent  (note  that  completely  different  pseudo-random  numbers  sequences 
i e  generated  by  the  algorithm  on  the  two  test  machines).  The  results  of 
I'abie  2  also  suggest  that  the  performance  of  SICMA  is  very  satisfactory 
row  the  point  of  view  of  dependability  (only  2  incorrect  claims  on  the 
’Marge"  dynamic  range  machine  when  Ng^g  >  3  and  on  the  "small"  dynamic 
range  machine  when  NgUC  >  4)  and  robustness  (no  overflows  on  both 
machines) . 

Unfortunately,  given  the  state  of  the  art  on  mathematical  software 
! or  global  optimization,  it  has  not  been  possible  to  make  conclusive  com¬ 
pel  i sons  with  other  packages. 

Finally,  we  note  that  a  smaller  value  of  Ng^  gives  a  much  cheaper 
method  (less  function  evaluations)  at  the  expense  of  a  loss  in  effective- 
'greater  number  of  failures). 
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